# ------------------------------------------------------------------------------#
# Description: Create Table 9: Robustness of peer effects to controls and sample criteria
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, modelsummary, ggplot2)
options(scipen = 999)

# Load adoption data -----------------------------------------------------------
load("data/adoptions/adoptions.RData")

# Merge peer adoption buckets --------------------------------------------------
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

# Merge peer eligibility buckets -----------------------------------------------
load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

# Merge total peers ------------------------------------------------------------
load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

# Merge parcel data ------------------------------------------------------------
load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

# Merge census data ------------------------------------------------------------
load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale     := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale  := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale      := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale      := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale    := peer.adopt.both.100 / 100]

# Keep post-2010 only 
dat <- dat[year > 2010]

# Create average eligibility variable ------------------------------------------
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Merge robustness datasets -----------------------------------------------------
load("data/shoreline/shoreline.RData")
setnames(near.shore, "near_dist", "near_shore")
dat <- merge(dat, near.shore[, .(pin, near_shore)], by = "pin", all.x = TRUE)
rm(near.shore)

load("data/arterial/arterial.RData")
setnames(near.arterial, "near_dist", "near_road")
dat <- merge(dat, near.arterial[, .(pin, near_road)], by = "pin", all.x = TRUE)
rm(near.arterial)

load("data/mandatory_gsi/peer_adopt_buckets_any_mandgsi.RData")
dat <- merge(dat, peer.adopt.bucket.any.mandgsi, by = c("pin", "year"))
rm(peer.adopt.bucket.any.mandgsi)

load("data/dem/all_parcels_dem.RData")
dat <- merge(dat, parcels.dem, by = "pin")
rm(parcels.dem)

# Create additional controls ---------------------------------------------------
dat[, close.road := ifelse(!is.na(near_road), 1, 0)]
dat[is.na(near_road), near_road := 0]

# Estimate models --------------------------------------------------------------

m1 <- feols(adopt.any*100 ~ peers.100 | year + blkgrp |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m2 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m3 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + near_shore |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m4 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob +
              close.road + close.road:near_road |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m5 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + peer.adopt.any.mandgsi.100 |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m6 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + elevation_ft |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m7 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + 
              near_shore + near_road + peer.adopt.any.mandgsi.100 + elevation_ft |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp", data = dat[elig.cs == 1 & post.rainwise == 0])

m8 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + 
              near_shore + near_road + peer.adopt.any.mandgsi.100 + elevation_ft |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0 & near_shore >= quantile(dat$near_shore, 0.2)])

m9 <- feols(adopt.any*100 ~ peers.100 + sqft + lot + beds + baths + yearbuilt +
              ren.pre90 + ren.90.99 + ren.00.09 + ren.10 + water.prob + 
              near_shore + near_road + peer.adopt.any.mandgsi.100 + elevation_ft |
              blkgrp + year |
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0 & near_road >= quantile(dat$near_road, 0.2)])

# Diagnostics and format extralines ---------------------------------------
load("data/iv/iv_diag_robust.RData")
for (i in 1:9) {
  assign(paste0('eff.f.', i), round(iv_diag_robust[[i]]$F_stat[4], 2))
  assign(paste0('ar.ci.', i), paste0('[', round(iv_diag_robust[[i]]$AR$ci[1], 2), ', ',
                                     round(iv_diag_robust[[i]]$AR$ci[2], 2), ']'))
}

# Variable names --------------------------------------------------------
varnames <- c("peer.adopt.any.100" = "$\\widehat{\\text{Peer Adoptions}}$")

# Preview table ---------------------------------------------------------------
etable(m1, m2, m3, m4, m5, m6, m7, m8, m9, keep = "%peer.adopt.any.100",
       fixef_sizes = FALSE, depvar = FALSE, se.below = TRUE, dict = varnames,
       fitstat = c("n", "ivfall"),
       extralines = list(
         "-Parcel Controls" = c("No", rep("Yes", 8)),
         "-Water Control"   = c("No", "No", "Yes", rep("No", 3), rep("Yes", 3)),
         "-Road Control"    = c(rep("No", 3), "Yes", rep("No", 2), rep("Yes", 3)),
         "-Mand. GSI Control" = c(rep("No", 4), "Yes", "No", rep("Yes", 3)),
         "-Elevation Control" = c(rep("No", 5), "Yes", rep("Yes", 3)),
         "-Exclude near Water" = c(rep("No", 7), "Yes", "No"),
         "-Exclude near Road"  = c(rep("No", 7), "No", "Yes"),
         "_F" = list(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5,
                     eff.f.6, eff.f.7, eff.f.8, eff.f.9),
         "_AR CI" = list(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5,
                         ar.ci.6, ar.ci.7, ar.ci.8, ar.ci.9)
       )
       )

# Save table -----------------------------------------------------------------
style_tex <- style.tex(
  tablefoot.value = '', model.title = '', depvar.title = '',
  fixef.title = '\\midrule \\emph{Fixed effects, controls, and sample}',
  var.title = '\\midrule', stats.title = '\\midrule')

etable(m1, m2, m3, m4, m5, m6, m7, m8, m9, tex = TRUE,
       keep = "%peer.adopt.any.100", fixef_sizes = FALSE, depvar = FALSE, se.below = TRUE,
       dict = varnames, digits = 2, fitstat = c("n", "ivfall"),
       style.tex = style_tex,
       extralines = list(
         "-Parcel Controls" = c("No", rep("Yes", 8)),
         "-Water Control"   = c("No", "No", "Yes", rep("No", 3), rep("Yes", 3)),
         "-Road Control"    = c(rep("No", 3), "Yes", rep("No", 2), rep("Yes", 3)),
         "-Mand. GSI Control" = c(rep("No", 4), "Yes", "No", rep("Yes", 3)),
         "-Elevation Control" = c(rep("No", 5), "Yes", rep("Yes", 3)),
         "-Exclude near Water" = c(rep("No", 7), "Yes", "No"),
         "-Exclude near Road"  = c(rep("No", 7), "No", "Yes"),
         "_F" = list(eff.f.1, eff.f.2, eff.f.3, eff.f.4, eff.f.5,
                     eff.f.6, eff.f.7, eff.f.8, eff.f.9),
         "_AR CI" = list(ar.ci.1, ar.ci.2, ar.ci.3, ar.ci.4, ar.ci.5,
                         ar.ci.6, ar.ci.7, ar.ci.8, ar.ci.9)),
       file = "output/tables/table_9.tex", replace = TRUE)
